1 Poisson Events

1.1 Libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))

1.2 Parameters and risk

censoredProb <- 0.000
timeSpan <- 10
timeInterval = 0.01
InitialPopulatoin <- 1000
ContBetaRate_1 <- 0.00005
BinBetaRate_1 <- 0.0005
BinBetaRate_2 <- 0.0002
BaselineHazard <- 0.00001
BinPrevalence1 <- 0.10
BinPrevalence2 <- 0.10

RContBetaRate_1 <- ContBetaRate_1
RBinBetaRate_1 <- BinBetaRate_1
RBinBetaRate_2 <- BinBetaRate_2
RBaselineHazard <- BaselineHazard

1.2.1 Simulate

#source("C:/Users/jtame/Documents/GitHub/RISKPLOTS/PlotRiskCategoriesAndTimetoEvent.R")

source("~/GitHub/RISKPLOTS/simulate_events.R")

1.3 RRplots Low Prevalence

plotTimeInterval <- timeSpan*3 #Not calibrated.

hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)


rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04534 0.11878 0.14200 0.24371 0.17643 0.95448

table(simulatedDataFrame$status)

0 1 897 103

par(cex=0.75)
RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=simulatedDataFrame$time,
                     title="Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=plotTimeInterval)

RRAnalysisCI_LP <- RRAnalysisCI

par(op)

1.3.1 By Risk Categories and Mean Time to Event

source("~/GitHub/RISKPLOTS/PlotRiskCategoriesAndTimetoEvent.R")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 103 84.1 124.9 107.0 0.962 0.786 1.17 0.735
low 51 38.0 67.1 49.5 1.031 0.768 1.36 0.776
90% 52 38.8 68.2 57.8 0.899 0.671 1.18 0.510 
pander::pander(timesdata)
  O:Low O:At Risk > 0.561 E:Low E:At Risk > 0.561
1Q 2.79 1.33 81.8 7.16
Median 4.94 3.68 100.5 9.66
3Q 7.04 5.95 123.9 14.39

1.3.2 Risk Calibration

op <- par(no.readonly = TRUE)


crdata <- cbind(simulatedDataFrame$status,pvalue,simulatedDataFrame$time)

#calprob <- CalibrationProbPoissonRisk(crdata,timeInterval=plotTimeInterval)
calprob <- CalibrationProbPoissonRisk(crdata)

( 8.214828 , 632.9083 , 0.9425186 , 103 , 106.8228 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.0275 0.29 8.21

1.3.3 After Calibration

h0 <- calprob$h0

caldata <- cbind(simulatedDataFrame$status,calprob$prob)

par(cex=0.75)

rrAnalysisTrain <- RRPlot(caldata,atRate=c(0.90),
                     timetoEvent=simulatedDataFrame$time,
                     title="Cal. Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

1.3.4 By Risk Categories and Mean Time to Event after calibration

source("~/GitHub/RISKPLOTS/PlotRiskCategoriesAndTimetoEvent.R")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 103 84.1 124.9 107.0 0.962 0.786 1.17 0.735
low 51 38.0 67.1 49.5 1.031 0.768 1.36 0.776
90% 52 38.8 68.2 57.8 0.899 0.671 1.18 0.510 
pander::pander(timesdata)
  O:Low O:At Risk > 0.561 E:Low E:At Risk > 0.561
1Q 2.79 1.33 81.8 7.16
Median 4.94 3.68 100.5 9.66
3Q 7.04 5.95 123.9 14.39

1.4 Second Simulation

ContBetaRate_1 <- RContBetaRate_1*3
BinBetaRate_1 <- RBinBetaRate_1*2
BinBetaRate_2 <- RBinBetaRate_2*2
BaselineHazard <- RBaselineHazard*3

source("~/GitHub/RISKPLOTS/simulate_events.R")

1.5 RRplots High prevalence

plotTimeInterval <- timeSpan # Calibrated

hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)


rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04534 0.11878 0.14200 0.22050 0.17605 0.88130

table(simulatedDataFrame$status)

0 1 782 218

par(cex=0.75)

RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=simulatedDataFrame$time,
                     title="Simulation 2",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=plotTimeInterval)

par(op)

1.5.1 By Risk Categories and Mean Time to Event

source("~/GitHub/RISKPLOTS/PlotRiskCategoriesAndTimetoEvent.R")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 218 190.0 249 216 1.008 0.879 1.15 0.892
low 99 80.5 121 109 0.911 0.740 1.11 0.388
90% 119 98.6 142 111 1.075 0.891 1.29 0.419 
pander::pander(timesdata)
  O:Low O:At Risk > 0.309 E:Low E:At Risk > 0.309
1Q 3.10 1.82 28.4 3.60
Median 5.33 4.21 33.1 5.11
3Q 7.41 6.27 37.7 8.15

1.5.2 Comparing the ROC

pander::pander(roc.test(RRAnalysisCI_LP$ROCAnalysis$ROC.analysis$roc.predictor,RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor))
DeLong’s test for two ROC curves: RRAnalysisCI_LP$ROCAnalysis$ROC.analysis$roc.predictor and RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor (continued below)
Test statistic df P value Alternative hypothesis AUC of roc1
-1.22 1698 0.222 two.sided 0.724
AUC of roc2
0.769

1.6 Third Simulation

ContBetaRate_1 <- RContBetaRate_1*10
BinBetaRate_1 <- RBinBetaRate_1*5
BinBetaRate_2 <- RBinBetaRate_2*5
BaselineHazard <- RBaselineHazard*5

source("~/GitHub/RISKPLOTS/simulate_events.R")

1.7 RRplots Higher prevalence

plotTimeInterval <- timeSpan # Calibrated

hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)


rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.09937 0.31029 0.36880 0.44629 0.44870 0.99573

table(simulatedDataFrame$status)

0 1 553 447

par(cex=0.75)


RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=simulatedDataFrame$time,
                     title="Simulation 3",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=plotTimeInterval)

par(op)

1.7.1 By Risk Categories and Mean Time to Event

source("~/GitHub/RISKPLOTS/PlotRiskCategoriesAndTimetoEvent.R")

pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 447 407 490 446 1.00 0.912 1.10 0.943
low 254 224 287 251 1.01 0.890 1.14 0.850
90% 193 167 222 195 0.99 0.856 1.14 0.943 
pander::pander(timesdata)
  O:Low O:At Risk > 0.450 E:Low E:At Risk > 0.450
1Q 2.12 1.35 10.1 1.59
Median 4.63 2.92 11.3 2.39
3Q 7.29 5.16 13.5 3.53

1.7.2 Comparing the ROC

pander::pander(roc.test(RRAnalysisCI_LP$ROCAnalysis$ROC.analysis$roc.predictor,RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor))
DeLong’s test for two ROC curves: RRAnalysisCI_LP$ROCAnalysis$ROC.analysis$roc.predictor and RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor (continued below)
Test statistic df P value Alternative hypothesis AUC of roc1
0.0419 1510 0.967 two.sided 0.724
AUC of roc2
0.723